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We propose a new theoretical approach to ground and low-energy excited states of nuclei extending 
the nuclear mean-field theory. It consists of three steps: stochastic preparation of many Slater 
determinants, the parity and angular momentum projection, and diagonalization of the generalized 
eigenvalue problems. The Slater determinants are constructed in the three-dimensional Cartesian 
coordinate representation capable of describing arbitrary shape of nuclei. We examine feasibility 
and usefulness of the method by applying the method with the BKN interaction to light 4iV-nuclei, 
12 C, 16 0, and 20 Ne. We discuss difficulties of keeping linear independence for basis states projected 
on good parity and angular momentum and present a possible prescription. 
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I. INTRODUCTION 

One of the goals in the microscopic nuclear many-body 
theory is the ab-initio nuclear structure calculation start- 
ing with a fixed Hamiltonian. Indeed, recent Green's 
function Monte-Carlo calculation with a bare nucleon- 
nucleon (NN) interaction is a milestone in this direction 
[H [1 [3 Q However, these ab-initio calculations are 
still limited to nuclei with mass number less than around 
twelve and to the lowest energy state for each parity and 
angular momentum. The no-core shell model |f| utilizes 
unitary transformation to accommodate the short-range 
correlations. These introduce effective NN interactions 
(without phenomenological adjustments), however their 
applications are also limited to light nuclei near the closed 
configurations. Systematic description of ground and ex- 
cited states in nuclei in a wide mass region requires de- 
velopment of a new computational approach. 

One of difficulties of the nuclear many-body problem 
is due to the strong short-range correlation. This for- 
bids a naive mean-field approach using the bare NN in- 
teraction. In order to overcome this difficulty, effective 
NN interactions have been extensively developed in his- 
tory of the nuclear theory 0. Namely, the short-range 
behavior of two-body correlation is effectively renormal- 
ized in the interaction. Then, nuclear many-body wave 
functions should take account of long-range correlations 
only Most of microscopic nuclear structure models adopt 
the effective interactions; e. g., n uclear mean-field models 
HH G3 El, shell models [ll [3, cluster models [3, 
etc. The success of these models indicates that a variety 
of low-lying modes of excitation are governed by noth- 
ing but the long-range correlations. In the present study, 
we aim for developing a new method to treat the whole 



correlation of long range beyond the mean field, utilizing 
the effective interaction for the mean-field models. 

To test our theory, we will apply the method to light 
N = Z even-even nuclei. There are a variety of nu- 
clear models for light nuclei. Yet, the existing models 
are not satisfactory in some respects. The shell model 
nicely describes spectroscopic properties of light nuclei. 
However, since the model space is truncated to a spe- 
cific shell, states very different from the ground state, 
which require a wider space, cannot be described ade- 
quately. A classic example is the second J 71 = + state 
in 16 O. Although this is the lowest-lying excited state 
of this doubly closed-shell nucleus, the shell model fails. 
The state has been successfully described by the cluster 
model. The nuclear cluster models have provided a fruit- 
ful description of many light nuclei, especially for those 
states lying close to the threshold. The antisymmetrized 
molecular dynamics (AMD), which was first utilized for 
study of heavy- ion collision |15|. is an extension of the 
microscopic cluster model successful to describe shell- 
model-like states as well 0, EJ- However, the model 
space of the AMD is a superposition of small number 
of Slater determinants whose orbitals are restricted to a 
Gaussian form. These models, the shell model and the 
AMD, are applicable to either light nuclei or those near 
the closed configuration. In contrast, the nuclear mean- 
field theory has been successful to describe nuclei of a 
wide mass region with a few parameters associated with 
the (density-dependent) effective nuclear force. The the- 
ory provides a reasonable description with a single Slater 
determinant for total binding energy, radius, and defor- 
mation of ground states. However, the superposition of 
multiple Slater is often required. For example, when the 
mean-field solution violates symmetries of the original 
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Hamiltonian, one should superpose many Slater deter- 
minants to restore the symmetry ("projection") Q. The 
symmetry restoration is crucial for many properties of 
light nuclei. Recently, we have developed a method of 
the parity (angular momentum) projection before (after) 
variation, and have shown that the mean-field model is 
capable of describing some typical cluster structures in 
light nuclei |l8j . In this paper, we intend to further ex- 
tend the method to a kind of complete calculation of the 
long-range correlations, to obtain energies and wave func- 
tions for the ground and low-lying excited states starting 
from a nuclear mean-field Hamiltonian. 

Our method has some resemblance to the generator 
coordinate method (GCM) 0, d and the Mote Carlo 
shell model (MCSM) [IJ. In the GCM, the generator 
coordinate is adopted a priori, under a certain physi- 
cal intuition, to describe specific long-range correlations; 
e.g., quadrupole and octupole correlation. In most prac- 
tical calculations, the coordinate is limited to one dimen- 
sion. Our method stochastically take into account all 
the important correlations. In the MCSM, basis states 
are stochastically generated and selected, then the diag- 
onalization of the Hamiltonian is performed in the space 
spanned by those states. This concept is very similar to 
ours, but we use the mean-field-model Hamiltonian and 
our model space is much wider than that of the shell 
model, 

The paper is organized as follows. In Scc.[Hl we present 
the outline of our method. Its details are described in 
the following sections; Selection of Slater determinants 
and the parity and angular momentum projection are 
shown in Sees. IIIII and IIVI respectively. In Sec. [V] we 
discuss configuration mixing calculation and how to avoid 
numerical instability caused by the overcompleteness of 
a selected basis and numerical errors. Then, we test the 
accuracy of our approach by taking 16 O as an example. 
In Sec. lVIl we compare numerical results of 12 C and 20 Ne 
to experimental data. The summary is given in Sec. I VI il 

II. FORMALISM 

In this section, we present the outline of our method 
to illustrate its essence. Roughly speaking, our method 
consists of three steps; 

(1) Generation and selection of Slater determinants (S- 
det's) important for ground and low- lying excited 
states. 

(2) Parity and angular momentum projection. 

(3) Configuration mixing (diagonalization of the 
Hamiltonian) . 

Each of these steps is not as straightforward as it first 
looks. For the step (1), since the short-range correlation 
is renormalized in the effective interaction, we should be 
careful not to adopt a S-det involving components with 
very high momentum. For the step (2), since there is no 



symmetry restriction on the wave function, we need to 
carry out the projection with respect to the full three- 
dimensional Euler angles. The diagonalization in (3) is 
cursed by the well-known overcompleteness problem of 
non-orthogonal basis and also by accumulated numerical 
errors in the step (2). We will present prescriptions to 
overcome these problems in Sees. IIIl| IIVI andlVl respec- 
tively. 

Following the prescription given in Sec. IIIII many S- 
det's are stochastically generated, then those important 
for low-energy excitations are selected. Single-particle 
wave functions in each S-det are represented in the three- 
dimensional (3D) Cartesian coordinate space without any 
symmetry restriction. The selected S-det 's form a basis 
set, {|<i> n ); n = lj • • • , N}. We then make parity and an- 
gular momentum projection for each S-det. The Hamil- 
tonian and the norm kernels for the fixed parity (jr = ±) 
and angular momentum (J, if), where K indicates its 
body- fixed component, are given by 

\ H M K '] = ^A H l ) P±P KK^n')- (1) 

Here, H is the many-body Hamiltonian with effective in- 
teraction, P^ is the parity projection operator, and P^k 1 
is the angular momentum projection operator. Finally, 
we solve the following generalized eigenvalue problem, 

E tfi* " EH±)N SlK<)9n>K> = 0. (2) 
n'K' 

If the space spanned by the set of the S-det's, {|$„); n = 
1, • • • , N}, is approximately complete for the long-range 
correlation, we should obtain a convergent solution for 
the ground and the low-lying excited states. Note that 
this is merely an outline of the method. As a matter of 
fact, in order to avoid the zero eigenvalues of the norm 
kernels, we will screen the selected basis states and mod- 
ify Eqs. This prescription will be given in Sec. IV! 

To demonstrate applicability of our method, we show 
numerical calculations employing the simplified mean- 
field Hamiltonian, the so-called BKN force |2l). The 
BKN force consists of two-body plus three-body forces. 
The two-body force consists of those of zero-range part 
(to), finite-range Yukawa part (Vy), and the Coulomb 
part (Vc). The three-body term is a zero-range force 
(t^). The Hamiltonian is given by 

H = E ++\ E y {2) m + 1 E ^ (3) m), 

i—l ^ ' ij ijk 

(3) 

where 

V {2 \ij) = t S{fi - rj) + V Y (n - f 3 ) + V c (n - f 3 ), 
V^{ijk) = t 3 6(n - f^Sifi - f fc ). 

The BKN interaction assumes that all nucleons have a 
charge e/2, and four nucleons with different spin and 



3 



--.(B) 




FIG. 1: Schematic picture of the energy surface. Two crosses 
represents minima of the energy surface. Three paths, (A), 
(B) and (C), show imaginary-time trajectories starting from 
different initial S-det's. The dotted arrow of (B) indicates the 
trajectory (B) passes through a shoulder state. 



isospin occupy the same spatial orbital. To express the 
orbital wave functions, we employ a grid representation 
discretizing 3D Cartesian coordinate. Grid points inside 
a sphere of 8 fm are adopted with grid spacing of 0.8 fm. 



III. PREPARATION OF BASIS 
A. Generation and selection of S-det's 

The first step is to prepare a set of S-det's that span 
the space necessary for low-lying states with the full long- 
range correlations. We make use of the imaginary-time 
method starting from initial configurations which are 
generated stochastically. The imaginary-time method is 
often used to obtain self-consistent solutions. Instead, we 
utilize it for generating many kinds of low-energy collec- 
tive surfaces. We pick up S-det's on the way to the self- 
consistent solutions before reaching minima, and employ 
them in the configuration mixing calculation. 

Figure ^ shows a schematic picture of the imaginary- 
time calculations starting from different initial configura- 
tions. The imaginary-time iteration has a property suit- 
able for generating the basis to calculate the long-range 
correlations. It quickly removes high-energy components 
of the wave function in a early stage of the iteration. The 
S-det is expected to rapidly fall onto a potential energy 
surface important for low-energy modes of excitation. 
This is the very property we want, because we should ex- 
clude S-det's which take account of the short-range corre- 
lation in the Hamiltonian. Therefore, we simply dispose 
all the S-det's generated in the first few hundred itera- 
tions of the imaginary-time evolution, then select S-det's 
after the rate of energy decrease becomes relatively slow. 

A series of S-det's generated with the imaginary-time 
calculation starting from an arbitrary initial state even- 



tually converge to a self-consistent solution; either the 
Hartree-Fock ground state (path (A) and (B) in Fig. 
or local minima solutions (path (C) in Fig. QJ. During 
the iterations, it sometimes happens that the configu- 
ration changes very slowly and the state stays almost 
unchanged for a long period of the iterations (a part pre- 
sented by the dotted arrow of path (B)). This is called 
a shoulder state. Although these shoulder states are not 
self-consistent solutions, they may play an important role 
for the low-lying excitation spectra and the ground-state 
correlation. 

We repeat the imaginary-time iteration many times 
starting from different initial configurations. We con- 
struct the initial S-det's by a stochastic procedure: The 
single-particle orbitals of the initial S-det are in a Gaus- 
sian form whose centers are randomly chosen. After gen- 
crating large number of imaginary-time trajectories, we 
may expect that those S-det's span the complete space 
for calculation of the long-range correlations. 

Figure [5] is an example of the actual imaginary-time 
calculations for 16 O, showing the energy expectation 
value, E(N it ) = ($(AT it )|iJ|$(AT it )), as a function of 
the iteration number, Nt- The path is similar to (B) 
in Fig. passing through a shoulder state. In iVjt < 
100, the energy expectation value decreases very rapidly. 
From Nn = 200 to 1500, the energy decreases very slowly, 
corresponding to a shoulder state. We have found that 
this shoulder state corresponds to the cluster structure 
of 12 C+a, which is considered as a dominant component 
of the first excited state of 16 O in the cluster model stud- 
ies. The dashed and dash-dotted curves are the energy 
expectation value after parity projection, E^(Nt) = 

($( ± )(7V it )|iJ|$(±)(iV it ))/(<f( ± )(7Vit)|<i>W(A r it)), where 
\& ± \N it ))=P ± mNit)). 



B. Selection of S-det's 

During the imaginary-time iterations of JVtotai steps, 
S-det's for every N s iterations are taken as candi- 
dates of the basis states. Thus, the S-det's at Nf t — 
N s , 2N S , • • • , k n N s are nominated first. The number of 
S-det's taken from a single path is k n = Noted/ N- The 
typical numbers are N s — 50 and iVtatal = 2000, leading 
to k n = 40. However, we cannot include all these S-det's 
in the basis set of the configuration mixing calculation, 
because too many S-det's lead to a numerical instabil- 
ity caused by the overcompleteness. Thus, we need to 
reduce their number. Here, we impose two additional 
constraints on those candidates: 



(a) E(Ng) < £ HF + 30 MeV - 

(b) Overlap between any pair of selected S-det's must 
be less than 0.7 (see below for details). 

The condition jsj means that the energy expectation 
value of each S-det, E(N? t ) = ($(A^ t )|tf |$(A^ t )), should 
not be so large because we are interested in low-lying 
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FIG. 2: An example of the imaginary-time evolution in 16 O 
started from a randomly generated S-det. Solid line indicates 
energy expectation value of the S-det, \&(Nit)), as a function 
of iteration number, Nit- Snapshots of the density distribu- 
tion at every 500 iterations are shown. The dashed and the 
dash-dotted line indicate energy of even and odd parity com- 
ponent of the S-det, respectively. The imaginary-time step of 
At = 0.001 h/MeV is adopted in the calculation. 



where R l ,i — 1 • ■ • 24 are special rotations and inversions 
corresponding to permutation of the axes (x,y,z). 

In practice, the Hartree-Fock state |$hf) is always 
selected first. Then, we start the examination of con- 
straints (jsj and © for generated S-det's in the 
ascending order of the energy expectation value. Here, 
we employ not the energy expectation values with re- 
spect to the S-det, but those with respect to the states 
projected onto negative parity, P~\^(N^ t )). This makes 
it easier to select exotic deformations which often ap- 
pear in the negative-parity excited states. For the case 
in Fig. states at w 2000 are examined first. Since 
those around Nu as 500 also show local minima in the en- 
ergy surface of negative parity, they are also nominated 
with high priority. 

The number of S-det's satisfying the criteria (jsj) and 
(|b| are typically from zero to three in a single trajectory 
of the imaginary-time evolution. Apparently, larger the 
number of selected S-det's, M, more difficult it becomes 
to find the (M + l)-th S-det to satisfy the condition (JbJ. 
In actual calculations, about 100 imaginary-time trajec- 
tories will be repeatedly generated to obtain about fifty 
S-det's which satisfy these conditions. 



states and the long-range correlations only. In the 
present work, we adopt the cut-off energy as 30 MeV 
above the Hartree-Fock ground-state energy. 

The second condition is directly related to the linear 
independence among the S-det's. In order to determine 
whether to include a candidate in the basis set, we exam- 
ine the overlaps between the new S-det (candidate) and 
all the S-det's which have been already included in the 
basis set. If the maximum of the absolute values of the 
overlap is less than a certain value, we add the candidate 
to the set of basis states. Since we make parity and an- 
gular momentum projection later, it is desirable to check 
this condition for projected wave functions. However, it 
costs too much to achieve angular momentum projection 
at this stage. Instead, we examine the overlap for differ- 
ent configurations produced by the interchange and the 
inversion of the Cartesian axes. These transformations 
correspond to 24 choices of the coordinate system, and 
are easily realized in the 3D Cartesian coordinate rep- 
resentation. The condition JEJ for adding a new S-det 
|$(A^)) to the selected basis set {!$„); n = 1, • • • ,M} is 
expressed by 

|($(JVa)|P ± fl < |$n)| 

K$(AS)|p±|$(JvS))($ n |p±|* B )|4 ' ' 

for n = 1, • • • , M, (4) 



IV. PARITY AND ANGULAR MOMENTUM 
PROJECTION 

A. Projection with respect to 3D Euler angles 

For each S-det in a set {|<6 n );n = 1, 2, • ■ ■ , N } which 
are generated and selected in Sec. IIIII we perform the 
parity and angular momentum projection. Since we con- 
struct wave functions employing the 3D Cartesian coor- 
dinate representation without any restrictions on the spa- 
tial symmetries, the full three-dimensional angular mo- 
mentum projection is necessary. The three Euler angles, 
(a, /3, 7), are abbreviated as f2. The angular momentum 
projection operator P^ K is given as usual by 

P J mk = J cMD J M * K (n)R(n), (5) 

where D'[ IK (fl) is the Wigner's D-function and R(fl) is 
the rotation operator. 

Quantities necessary for solving a generalized eigen- 
value equation in Sec.^are matrix elements of the norm 
and Hamiltonian kernels. For two S-det's, |$„) and |$ n '), 
we calculate matrix elements for parity and angular mo- 
mentum projected wave functions. 
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nK,n'K' 
H J(±) 
n nK,n'K' 



2J+ 1 



]( _ , <MD J K * K ,{n) x ($ n |e" iQ J * ^ (l±P)e-^e-^|$ n ,), 



(6) 



where K and .K 7 specifies quantum numbers for body- 
fixed component of the angular momentum operator. 
Since the rotation operator commutes with the Hamil- 
tonian, we apply rotations of angles [3 and 7 to the ket 
vector |<£v) and the rotation of a to the bra vector ($„| 
in the practical calculations. The parity projection oper- 
ator is given by P = (1/2)(1±P), where P is the space 
inversion operator. 

The parity transformation and rotation of the S-det, 
|$) = det{|0fe)}/\/]4T, are achieved by the correspond- 
ing transformation of their single-particle orbitals, \4>k), 
k = 1, • • • , A. The rotation of finite angle is realized by 
successive rotations of a small angle. For instance, for a 
rotation of angle a around z-axis, 

|^) ee e- ia % fc ) = (e- <Ao &) JVdlT \4> k ), Aa = a/N div . 

Each small-angle rotation is performed using the Taylor 
expansion of the rotation operator; 

« "ft^ie, <b) 

fc=0 

where A^ max = 4 gives an accurate result. We usually 
employ Aa = 27r/360. 

The integrand of Eq. JBJ is the overlap/Hamiltonian 

matrix element between two different S-det's, e~ taJ - 

and (p) e ~ J ' 3 ' / f e _l7l7z |$„'). These matrix elements are 
simply expressed in terms of the interstate density matrix 
defined in Appendix. 

B. Numerical details of the projection 

We here discuss numerical accuracy of the 3D angular 
momentum projection. Numerical error in the matrix 
elements of Eq. © may cause a serious trouble when 
we solve the generalized eigenvalue problem of Eq. J2J. 
For instance, the norm matrix, should be pos- 

itive definite, however, in practice, calculated norm ma- 
trix suffers from many negative eigenvalues though their 
absolute values are small. 

The finite difference approximation for the angular mo- 
mentum and the finite-order expansion for the rotation 
operator in Eq. JSJ are good approximation. We have ex- 
amined the identity of the single-particle orbitals before 



and after rotating over 2ir. The overlap between these 
two single-particle wave functions is very close to unity, 
with error less than 10~ 4 . Therefore, the error in each ro- 
tated wave function is relatively small. However, it seems 
that these numerical errors are accumulated during the 
3D integration over Euler angles, < a < 2tt, < (3 < tt, 
and < 7 < 2tt. 

The numerical integration is carried out using the 
trapezoidal rule with the finite-number discretization. 
Figure |3 shows the norm eigenvalues calculated for J 71 " = 
3~ in 16 0. Fifty Slater determinants are generated in the 
procedure explained in Sec. lIIII thus the dimension of the 
norm matrix is 350 (— 7 x 50 where 7 is the number of 
different K quantum numbers). The eigenvalues of norm 
matrix N n y K ' n , K , are plotted in descending order. The 
left three panels show the eigenvalues when the number 
of grid points of the angle /3 is varied from 20 to 80 while 
those for a and 7 are fixed at 20. For the right panels, we 
vary the number for a and 7 from 10 to 40, being fixed at 
20 for (3. The norm eigenvalues are converged for a and 
7 if we adopt 15 or more grid points for their discretiza- 
tion. In contrast, the convergence with respect to (3 is 
rather slow and we still have about 70 negative eigen- 
values with 80 grid points. Apparently, it is desirable to 
have larger number of grid points for the discretization 
of (3 angle. However, we must make a compromise and 
sacrifice some accuracy, because the calculation for the 
Hamiltonian kernels in Eq. © are very demanding in the 
3D rotation. In the present work, we employ discretiza- 
tion of 15 grid points for a and 7 and 20 points for (3. 
Even with this discretization, we need to evaluate 4,500 
matrix elements for each pair of the S-det's. 

A price of sacrificing the accuracy is a complication 
of solving the eigenvalue problem, due to appearance of 
negative norm. In Sec we will explain how to cope 
with this difficulty. 



V. CONFIGURATION MIXING; ENERGY 
SPECTRA OF 16 

A. Zero- and negative-norm problem 

Although we collect a set of linearly independent S- 
det's in the procedure explained in Sec. 11111 the linear 
independence is often lost after the parity and angu- 
lar momentum projection. This will lead to number of 
eigenvalues of the norm matrix close to zero. Moreover, 
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FIG. 3: Eigenvalues of norm matrix for J™ = 3~ in O. The horizontal line indicates the sequential number according to the 
magnitude of the eigenvalues. The top, middle and bottom panels show those of No. 1 to 200, 200 to 320, and 320 to 350, 
respectively. The norm matrix is calculated with various discretization on Euler angles. In the left panels, f3 is discretized into 
20 (squares), 40 (circles), and 80 (triangles) points, while a and 7 are discretized into 20 points. In the right panels, a and 7 
are discretized into 10 (squares), 20 (circles), and 40 (triangles) points, while /3 is into 20 points. 



the numerical error in the angular momentum projection 
even produces the negative eigenvalues. In order to solve 
the eigenvalue equation J5J), we must remove states which 
cause the zero and negative eigenvalues. In this section, 
we give a possible prescription for this. 

For each angular momentum state J, we reduce the 
dimension of the configuration space as follows. We first 
diagonalize the norm matrix in different K states for each 
S-det. This is the diagonalization of (2 J + 1) x (2 J + 1) 
matrix, 

J 

1S nK,nK' V K' ~ A nv V K • l y I 

K'--J 

Here the eigenstates are labeled by v. After the K- 



diagonalization, the basis is specified by the label in- 
stead of K. The collected basis states are denoted as 

M =J2 V K P MKP ± \*n), M = — J, • • • , J. 
K 

(10) 

The magnetic quantum number in the laboratory frame, 
M, is a trivial conserved quantity simply giving the 
(2 J+l)-fold degeneracy for each state, thus omitted here- 
after. At this stage, we exclude states whose eigenvalues 
Xni^ are smaller than 10~ 2 . A small eigenvalue will 
appear when, for example, the S-det is an approximate 
eigenstate of a specific K quantum number. In the case of 
J v = 3~ of 16 O, we have 350 configurations constructed 
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from fifty S-det's, among which about one hundred con- 
figurations have Xnu less than 10~ 2 and are discarded. 

Although the states in {\^n[^)}n,v,j are truncated ac- 
cording to the norm eigenvalues, in order to avoid numer- 
ical instability, we need to reduce the number further- 
more. For this purpose, we consider a following "normal- 
ized" norm matrix, 



corresponding to an eigenvalue Aj as 

I*/ m(±) )m = -±= £ -^k= £ v%p J MK p±\s n ). 

(13) 

Using these states as a basis, we calculate the norm and 
Hamiltonian kernel matrices, 



N 



J(±) 



j(±) 



nu.n' v' 



N J(±) 



1/2 



(V, (±) , , 



1/2' 



(11) 



and solve the generalized eigenvalue equation 



(14) 



which is constructed so as to make the diagonal elements 
Nninl equal to unity. We make a further selection ac- 
cording to the magnitude of eigenvalues of this matrix, 
Ai, obtained by solving 



ViV J(±) u 



Tl'V 



= \u 



'nvi 



(12) 



n; 



J(±) _ p j(± 



y J y j 



0. (15) 



We obtain the energy eigenvalues E J ^ and the eigen- 

, ~J(±) 
vectors ^ . 



for each parity and angular momentum sector. The exis- 
tence of small eigenvalue, A, -C 1, indicates a strong over- 
completeness of the basis set. We impose the condition 
on the eigenvalues, that the smallest eigenvalue, A m j n , 
must be greater than 10 -3 . This is done by the following 
procedure. First, we calculate eigenvalues of 2 x 2 matrix 

of ^£l< for a11 P° ssible P airs 01 (l*^), l*^))- H 
the smaller eigenvalue is less than 1CP 3 , we remove one of 
them according to the magnitude of its diagonal element 
(remove | ¥„,„,) if N^] n , vl < Nnlfnv)- The number of 

basis states, {l^n^)} surviving these screenings is now 
denoted as N sc . For the J 77 = 3~ states of 16 0, N sc is 
of order of one hundred. If the A m i n of the N sc x N sc 
matrix, N^^, v ,, is larger than 10 -3 . we can proceed 
to the configuration mixing calculation to solve Eq. 1)150. 
Otherwise, we will further reduce the number of states: 
We diagonalize the matrix N^J, V , in a space spanned by 

the basis except for a single state, \^fm^). This is the 
diagonalization of the (N sc — 1) x (N sc — 1) matrix. We 

do this N sc times for all possible \^mp), in order to find 



J(±)\ 

mLi /ex 



, for which the minimum eigenvalue 



the one, |* 

of the remaining (N ac — 1) x (N sc — 1) matrix becomes 

the largest. This state, |^jn/» )ex> is removed from the 
basis set. This process is repeated and the number of 
basis is reduced one by one, until all the eigenvalues, Ai 
(i = 1, • • • , -Wjf ), becomes larger than 10~ 3 . In the 
case of J 7 ' = 3~ in 16 O, several dozens of configurations 
are discarded in this screening, so that the 



of I* 



3(-)\ 
mil 



3(— ) 

final number of the basis states is N b « 50. Note 
that is the number of states in {|^nS )} 5 thus 

the number of adopted S-det's (that of |$ n )) is i n gen- 
eral less than iv/^. 

In order to check numerical accuracy and stability, it 
is convenient to define a "normalized" norm eigenstate 



B. Quality of solutions 

In this section, we examine quality of solutions ob- 
tained by diagonalizing Eq. I|15|) and how "complete" the 
selected basis is. Let us emphasize again that we do not 
intend to obtain the exact eigenstates of a given Hamilto- 
nian. The exact ground state of a Hamiltonian with the 
zero-range interaction such as Eq. @ perhaps leads to 
an unphysical solution. Instead, we aim to take into ac- 
count correlations of its long-range part only. Therefore, 
we examine whether the method can produce convergent 
results for low-lying states. 

Let us suppose that we have selected N^^' basis states 
for specific parity and angular momentum, (J, M, ±). We 
solve Eq. ifT^ to obtain eigenvalues, At, and vectors, u l nv , 
then construct the norm eigenstates of Eq. (|13fl . |^^ ± ' ) ). 
First, the states {I'J'f )} are sorted according to the 
magnitude of their norm eigenvalues^ A^. Thus, the eigen- 
states are arranged in sequence of Ai > A2 > • • • > Ajv b - 
The middle and bottom panels in Fig. 0] show distribu- 
tions of Ai and the diagonal elements of the Hamiltonian, 
E J % {±) = H-} ±] , respectively, for J 71 = 0+ (left) and 3" 
states (right) in 16 0. In this calculation, there are 19 basis 
states for J 77 = + and 53 for J 77 = 3~ for which all the 
norm eigenvalues are larger than 10~ 3 . The Ai decrease 
linearly in the logarithmic scale. An interesting thing is 
the fact that the energy expectation values, , are 

closely correlated with the norm eigenvalues Ai. The en- 
ergies roughly show a monotonic increase with i, 
as Ai decrease. This may justify the screening process 
to discard states with small norm eigenvalues, because 
those states possess large energy expectation values and 
are expected not to play a significant role for low-energy 
excitations. 

The top panels in Fig.0]show resultant energy eigenval- 
ues, E J ^\ obtained by the diagonalization of Eq. \lb\. 
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FIG. 4: Energies and norm eigenstates of J 7r =0 + states (left) and 3 states (right) in 16 0. The top panels show calculated 
energy E J ^ of Eq. I15H as a function of dimension of the adopted model space. The middle and bottom panels show Xi and 



E 



j(±) 



H J U , respectively. See text for details. 



The horizontal axis indicates the number of basis states 
included in the calculation, which is increased one by 
one from left to right; {|*f (±) )}, {|*i (±) >, |* 2 7(±) )}, ■ ■ ■ , 
{|*i (±) ), • • • , l*^}}. Calculated spectra for ,V = 0+ 
and 3~ states in 16 O are shown in the left and right 
panels, respectively. As is seen in the figure, low-energy 
spectra become almost invariant with respect to the in- 
clusion of new basis states. In other words, energies of 
the ground and low-lying states are insensitive to the 
inclusion of states with small norm eigenvalues. These 
convergent behaviors suggest that the long-range corre- 
lations for low-lying states are taken into account in the 
calculation. 

To further examine the completeness of the basis 
states, we check the identity of results produced with 
different sets of basis states (initially generated with dif- 
ferent random numbers). If our prescription provides a 
complete set of basis for the long-range correlations of 
the Hamiltonian, the energy spectra should not depend 
on the initial S-det's from which the imaginary-time it- 
eration started. In Fig. [SJ excitation energies in the 16 
nucleus calculated with four different sets of S-det's are 



compared. In these four independent calculations, differ- 
ent seeds for the random number were used in preparing 
the initial S-det's. The energies of the lowest and the 
next lowest states for each parity and angular momen- 
tum (J*) coincide to each other in reasonable accuracy. 
For example, negative-parity excitations of l - , 2~ , and 
3~ states appear below 15 MeV, and there are no other 
states in this energy region. The results become less re- 
liable for higher states in each J 7r sector. The second + 
state (first excited + ) appears around 15 ~ 17 MeV in 
all calculations. However, excitation energy of the third 
(second excited) + state in the bottom-left panel is no- 
tably higher than the other three. 

The excitation energy of the second + state is much 
higher than the experimental value (6.05 MeV). The 
three negative-parity excited states with J 71 " = 3~, 1 _ , 
and 2~ , are observed at 6.13, 7.12, and 8.87 MeV, respec- 
tively. Since the BKN interaction adopted in the present 
work, which does not contain the spin-orbit interaction, 
is too simple to give a quantitative description of nuclear 
structure, one should not seriously take these discrepan- 
cies between the calculation and the experiment. 
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FIG. 5: Excitation energies of le O. Symbols at the bottom of each panel indicate the quantum numbers, J w . Results of 
calculations employing four different sets of Slater determinants are displayed. 



We also perform the same examination for other nuclei 
discussed below, 12 C and 20 Ne. The final number of basis 
states JV^ and behavior of the convergence is similar 
to those of 16 0. The comparison of results among sets of 
basis generated with different random numbers may pro- 
vide a useful information about reliability of calculations. 
From these analysis, we may judge how many eigenstates 
in each J T sector can be trusted. 

Before finishing this section, let us comment on the cut- 
off value of the norm eigenvalues. In the bottom panels of 
Fig. ^] the energy expectation values become somewhat 
scattered as the norm eigenvalues, Ai, approach to 10~ 3 . 
This may be an indication of the numerical instability. 
In fact, if we include basis functions with Xi < 10~ 3 , 
the configuration mixing leads to unphysical solutions. 
For instance, if we set a cut-off of A; > 10~ 4 , the ground 
state becomes completely different from the Hartree-Fock 
state and its energy is unreasonably lowered by the diag- 
onalization of the Hamiltonian. It is most likely that this 
problem originates from the numerical error in evaluating 
the matrix elements, especially in the angular momentum 
projection. 



VI. ENERGY SPECTRA OF 12 C AND 20 Ne 

The BKN interaction of Eq. (j3Jl is adopted for test- 
ing our new method. As we have mentioned before, one 
should not expect a quantitative description of the low- 
lying spectra. However, the present calculation gives a 
reasonable description for some excited states of light 
nuclei, especially for those composed of the LS'-closed 
clusters. In this section, we present calculated energy 
spectra of 20 Ne and 12 C nuclei. In these nuclei, there 
appear the LS'-closed clusters in the ground and excited 
states (a+ 16 for 20 Ne and 3a for 12 C). In this paper, 
we restrict our discussion on the energy spectra in these 
N = Z even-even nuclei. A detailed discussion on the 
structure of excited states including information on the 
transition matrix elements will be given in our next work 
employing a realistic Skyrme interaction. 

Figure El shows calculated energy spectra (left panel) 
and those in measurement (right) in 20 Ne. It has been 
well-known that the K v — + ground state band and the 
A' 71 " = 0~ negative-parity band starting with l - state at 
5.785 MeV constitute a kind of inversion doublet band 
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FIG. 6: Calculated and experimental excitation energy spec- 
tra of 20 Ne. Symbols at the bottom of each panel indicate the 
quantum numbers, J". 



FIG. 7: The same as Fig. ©but for 12 C. 



VII. SUMMARY 



of the 



cluster structure. This inversion doublet 



bands are reasonably described in our calculation. In 
the measurement, the lowest negative parity band is the 
K 7 * = 2~ band at 4.968 MeV. In our calculation, it is 
around 9 MeV, reflecting the importance of the spin-orbit 
splitting of p and d orbitals for this excitation. 

We next discuss results for 12 C whose spectra are 
shown in Fig. [7| Again, the calculation (left panel) is 
compared with measured spectra (right). The calcula- 
tion produce the ground state rotational band, but its 
moment of inertia is significantly larger than the observed 
values. In the negative parity, our calculation produces 
3~ state in the lowest energy, also 1~ and 2~ states at low 
excitation energies. These are qualitatively in agreement 
with experiments. In the positive parity excited states, 
the calculation indicates two + states around 10 and 
12 MeV. These may correspond to the measured states 
around 8 and 10 MeV. 

In Refs. [2^ , we have reported the variation after par- 
ity projection calculation employing the BKN interac- 
tion. The angular momentum projection after variation 
is achieved in the calculation employing Skyrme force 
[IM |24| . A part of results presented in this section 
coincide fairly well with these variation after parity pro- 
jection calculations. This suggests that the variation af- 
ter parity projection gives a dominant correlation for a 
certain class of states (with clustering) in 12 C and 20 Ne. 



In this paper, we report our attempt to develop a new 
computational method to include all the long-range cor- 
relations beyond the mean-field approximation. We aim 
at a systematic description of the ground and low-lying 
excited states using a mean-field Hamiltonian, without 
assuming their structure a priori. 

First, we stochastically generate many Slater determi- 
nants. The single-particle orbitals are expressed on the 
three-dimensional Cartesian grid representation. In or- 
der to remove high-energy components in those states, 
we use the imaginary-time iteration method. The 
imaginary-time evolution produces many trajectories im- 
portant for low-energy modes of excitation. We select 
some of these states to keep the linear independence. We 
then project them on good parity and angular momen- 
tum, and perform configuration mixing calculation. The 
BKN interaction is utilized to examine feasibility and 
difficulty of the method. We have found that there is 
a numerical difficulty to achieve the configuration mix- 
ing calculation. The eigenvalues of the norm matrix can 
be close to zero, when the selected states are overcom- 
plete. In the practical calculations, a small numerical 
error in the angular momentum projection results in the 
occurrence of negative eigenvalues of the norm matrix. 
We eliminate states responsible for these zero and nega- 
tive eigenvalues before solving the generalized eigenvalue 
problem. 
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We show calculated results for some light 4Y nuclei, 
12 C, 16 0, and 20 Ne. In these nuclei, appearances of vari- 
ous cluster states are known in excited states. Our calcu- 
lation provides reasonable excitation energies for a+ 16 
states of 20 Ne and 3a states of 12 C for which the spin- 
orbit interaction, which is not included in the BKN force, 
does not play an important role. 

The results calculated with different sets of random 
numbers are approximately identical to each other. How- 
ever, the discrepancy becomes more evident for states at 
higher energies. In the present level of accuracy, we may 
predict the lowest and possibly the second lowest states in 
each parity and angular momentum sector. Improvement 
in numerical accuracy, especially in the three-dimensional 
angular momentum projection, is desired for future work 
with more realistic interactions. 
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APPENDIX: MATRIX ELEMENTS BETWEEN 
NON-ORTHOGONAL SLATER DETERMINANTS 

In this appendix, we present useful expressions for cal- 
culating matrix elements such as the integrand of Eq. . 
In general, we discuss transition amplitude of an op- 
erator 6 between two different S-det's, ($|0|*). The 
S-det's are expressed in terms of orthonormal single- 
particle orbitals, |$) = det{\4>i(j))} / y/Al and |*) = 



det{|^i(j))}/V3!, with {(j>i\(f>j) = (ipi\i>j) = S i:j . Here, 
we assume |<fr) and \*Sf) are not orthogonal to each other. 
To calculate these matrix elements, it is convenient to 
define following orbitals, 

3 

where the matrix B is defined by B^ = (4>i\ijjj). The 
overlap matrix element is given by the determinant of B, 
($|^) = det-B. It can be easily confirmed that are 
bi-orthogonal to \4>j), having 

{<f>S j )=S ij . _ (A.2) 

We also note that the S-det constructed from \ipi) is pro- 
portional to \^), 

Since the S-det \^>) is now represented by single-particle 
orbitals which have a bi-orthogonal property (|A.2fl . 
the matrix elements ( < I > |0| x I') can be expressed in a fa- 
miliar form very similar to that of the expectation value, 
($|0|<f>). Suppose that, the expectation value of the op- 
erator O, ($|0|<i>), is expressed as a functional of the 
density matrix, Pij(<&) = ($|c]ci|$), 

0\p(*)] = (*|0|*). (A.4) 

Then, the matrix element between |$) and \^f) is ex- 
pressed as 

($|6|#) = ($|0|#)($|*) = 0[ j 5($*)]($|*), (A.5) 
where the interstate density /5(<I>^) is defined by 

= ($|ctci|§>. (A.6) 

Therefore, the matrix element between two S-dct's, 
($101^), is given by 0[p] x dct B, where the density ma- 
trix p is replaced by the interstate density matrix (|A.6(I . 
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